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Abstract- We propose a new method to numerically calculate higher-order correlation 
functions of primordial fluctuations generated from any early-universe scenario. Our key- 
starting point is the realization that the tree-level In-In formalism is intrinsically separable. 
This enables us to use modal techniques to efficiently calculate and represent non-Gaussian 
shapes in a separable form well suited to data analysis. We prove the feasibility and 
the accuracy of our method by applying it to simple single-field inflationary models in 
which analytical results are available, and we perform non-trivial consistency checks like 
the verification of the single field consistency relation. We also point out that the ie 
prescription is automatically taken into account in our method, preventing the need for 
ad-hoc tricks to implement it numerically. 
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1 Introduction 

The deviation from perfect Gaussian statistics of the primordial cosmological fluctuations 
enables us to discriminate amongst the candidate physical mechanisms that produced 
the seed primordial fluctuations (for recent reviews, see for instance [1-5]). With the 
imminent arrival of Planck data together with those from large scale structures surveys, 
it is fundamental that theorists have firm predictions in wide classes of models to match 
against these increasingly precise data. However, most non-Gaussian shapes up to now 
have been determined under restrictive approximations, typically of the slow- varying type, 
to enable an analytical calculation. Interesting inflationary dynamics nonetheless exist 
which violate these assumptions and for which no theoretical prediction is yet available, 
for instance models with a strong background time-dependence or multifield models with 
non-trivial trajectories (see for example [6] for a recent illustration). 

Procedures to numerically calculate the non-Gaussian properties of early-universe 
models do exist [7-10], but they rely on direct fc-space configuration by fc-space configu- 
ration calculations requiring tricks to implement the ie prescription necessary to project 
onto the interacting vacuum, like the introduction of a 'damping' factor into the required 
time-integral to ensure early-time convergence. 

In this paper, we propose an alternative method based on the separability of the tree- 
level In-In formalism and the modal decomposition techniques of Fergusson, Shellard and 
collaborators [11-20] (see also [21-23]). It enables one to numerically compute tree-level 



higher-order correlation functions of cosmological observables from any early-universe sce- 
nario, with three main advantages: 

• no ad-hoc trick is required to implement the ie prescription, which is automatically 
taken into account in our method. 

• the results of our procedure are smooth functions of their arguments, the fcj's, and 
not values of the relevant correlation functions on a grid of points in Fourier space. 

• the resulting shapes are directly represented in a separable form, and hence can 
be efficiently constrained by data, for instance from the Planck satellite. 

The outline of the paper is as follows: in section 2, we present the key-idea behind 
our method, namely the use of the intrinsic separability of the tree-level In-In formalism, 
and details about its numerical implementation. In section 3, we demonstrate its efficiency 
by calculating the tree-level bispectra generated by various operators in single-field test 
scenarios, like chaotic and Dirac-Born-Infeld inflation, and comparing them to known 
analytical results. We also perform non-trivial checks of our formalism, including the 
verification of the single field consistency relation between the scalar spectral index and 
the squeezed limit of the primordial bispectrum [24, 25]. We state our conclusions in 
section 4. 



2 Method 



2.1 Separability of the In-In formalism 

The central starting point of our method relies on the realization that the In-In formalism 
used to calculate primordial cosmological correlation functions [26] is intrinsically separa- 
ble, at least at tree-level. 

The In-In (also known as Keldysh-Schwinger) formalism gives the following prescrip- 
tion for calculating the expectation value of any operator Q(t) in the interacting vacuum 
of the theory: 



m)) = <oi 
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where |0) denotes the vacuum of the free theory, T the time-ordering product, Hj the 
interacting Hamiltonian, Q T (t) is in the interacting picture, and the ie indicates that the 
time integration contour should be slightly rotated into the imaginary plane in order to 
project onto the interacting vacuum [24]. 

Concentrating, for simplicity of presentation, on the bispectrum of the curvature 
perturbation ( generated in single-field models of inflation, one is led to evaluate, at tree- 
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level: 

(Ck^Kk^rKkatT)) = -i f dr'a(r')(0|C kl (r)C k2 (r)C k3 (r)F ( 3)(r')|0) + c.c. (2.1) 

J — oo(l— ie) 

where r denotes the conformal time such that dr = dt/a(t), H^) is the cubic interacting 
Hamiltonian, and all operators on the right-hand side are in the interacting picture (we 
omit the I for brevity). In practice, the dimensionless quantity we will compute and which 
is ultimately the most useful for data analysis is the so-called shape S(ki, fc 2 , £3), defined 
such that [12, 27] 

(C kl C k2 CO - (2vr) 7 5(k 1 + k 2 + k 3 ) S f^^ Al , (2.2) 

where A s is a normalization factor related to the amplitude of scalar fluctuations. Il- 
lustrating this on the two cubic interactions D J dtd 3 xa(t)g(t)((d() 2 and D 
f dtd 3 xa(t)g(t)((d() 2 , where g(t) denotes a generic time-dependent coupling constant, 
we find respectively that 

S(( dc) 2(ku k 2 , k s ) oc ik 2 • k 3 / dr' a{r')g(r') klCkA^Ck^' 

J-oc(l-ie) \ 

k 2 2Ck 2 (r)C k2 (r')) I k!Ck s (r)C k3 (r') J + c.c. + 2perms. (2.3) 



X 



and 



S m)2 (h,k 2 , k 3 ) oc *k 2 • k 3 [ T dr'a 2 (T')g(r') [ kl( kl (r)C kl (r') 

J-oo(l-ie) \ y 

x [klCkAr)C k2 (r') j ^3 2 a 3 (r)a 3 (^) j + c.c. + 2 perms. (2.4) 

where (' k denotes the time derivative of the mode function with respect to r, and k 2 • k 3 = 
\ {k\ — k\ — fcf) using the triangle condition ki + k 2 + k 3 = 0. 

The key-point to notice is that, up to explicitly known /c-dependent factors, like k 2 • k 3 
in Eqs. (2.3)- (2.4), the primordial shape generated by any cubic operator is explicitly given 
in a separable form, i.e. as the time integral of a product of three functions of respectively 
fcij ^3- This crucial fact therefore enables one to apply efficient separable techniques to 
determine a modal decomposition of the non-Gaussian profile S generated by any cubic 
operator. 
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2.2 Modal decomposition 

Following the modal techniques of Fergusson, Shellard et al, we wish to represent any 
primordial shape S(ki, fc 2 , k 3 ) (or similar quantities related to the relevant shape by mul- 
tiplication by known fc-dependent factors and/or symmetrizations) as 

S(k u k 2 ,k 3 ) = ^2a n Q n (k 1) k 2) k 3 ) (2.5) 

n 

where the Q n 's are suitable basis functions that are separable and orthonormal, so that 
the a n 's can be computed efficiently. For this purpose, we should halt and define both 
the space over which we want to expand the shapes and the scalar product with respect 
to which orthonormality is defined. Physically, the shapes S are defined on a domain in 
Fourier space such that: (i) all fc^'s G [fc m inj ^max], where fc min and /c max stand for respectively 
the largest and smallest scales of interest (which depend on the particular probe we are 
interested in, for example CMB or large scale structures), (ii) the wavenumbers are such 
that the triangle condition ki + k 2 + k 3 = holds, namely 

k\ < k 2 + k 3 for k\ > k 2 , k 3 , (and permutations) . (2.6) 

However, it is easy to realize that calculating the coefficients a n 's on such a tetrahedral 
domain would come at the cost of losing separability. Therefore, we consider instead the 
expressions (2.3)- (2.4) (and similarly for other vertices) as mathematically defining shapes 
over the cubic domain C = [k m [ ni fc max ] 3 , irrespective of the triangle condition. Of course, 
for physical interpretations, the shapes we will thus construct should be restricted in the 
end to the tetrahedral domain, but the use of shapes defined on C enables us to use the 
intrinsic separability of the In-In formalism. Indeed, using the canonical scalar product 
on C: 

F{S, S') = J S(k u k 2) h)S\k u fc 2 , fc 3 ) dfci dk 2 dk 3 , (2.7) 

orthonormal and complete families of separable basis functions Q n with the correct symme- 
try properties can be constructed out of products of Legendre polynomials. For instance, 
a completely symmetric function of (£4, fc 2 , £3), like the building block of Eq. (2.4): 

S m) <h,k 2 ,k 3 ) = [ T &T'a\T')g{T')(kl(; kl {T)e kl {T') 

J-oo(l-ie) \ 

x f klC k2 {r)C k2 {r') ) f k 2 3 C k3 (r)C k3 (r') ) , (2. 
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can be represented on C as in Eq. (2.5), using the basis functions 

_ nL(2^ + i) i / 2 A 

( U _ h. \3/2 ^nwns 



2k\ (fc m in + ^max)\ ^ ( 2&2 (^min + ^max) \ ^ / 2&3 (^min + ^max) 




I _L \- "111111 ' '"llld^/ 1 y-) 

x p («i I — i — tt" J p ^ \ — i — znr. ) \ k -k 

where Pi denotes the z-th order Legendre polynomial, the index n represents the symmetric 
triplet {ni, n 2 , n 3 }, and where A nin2n3 equals 1, v^, v6 if all the r^'s are equal, two of them 
are equal or all of them are different respectively. Using these orthonormal and separable 
basis functions, together with the separable representations of shapes as given by the In-In 
formalism, the modal decomposition coefficients a n 's in Eq. (2.5) can be computed as 
a n = f c SQ n , where the triple integral over the cube C can be crucially decomposed as 
products of three one-dimensional integrals over fci, k 2 and fc 3 respectively: 

/T I ^max 

dr'a 2 (r')g(r')[ / P ni ^WC(^) 
-oo(l— ie) \ J Li T1 

^ P n2 k 2 2 Ck 2 (r)C k2 (r') dk 2 U I P n3 klCk 3 (r)C k3 (r') dk 3 ) , (2.9) 

J \J ^min J 

and where the argument of the P n 's are (2h n — (fc m in + fc ma x))/(fcmax — hmm)- Of course, a 
similar expansion holds for non-symmetric functions, like the building block of Eq. (2.3) 

S t{d<)2 (k u k 2 , k 3 ) = f dr'a(r')g(r') f klC kl {r)C k [{r') 

J — oo(l— ie) 

x ^ 2 2 a 2 (r)C 2 (r')J [k!( k3 (T)Q 3 (T') j , (2.10) 

which is symmetric only in (k 2) ks)- The only difference in such cases is that the basis 
polynomials should only be symmetric in (k 2 ,ks). Following these methods, the shape 
generated by any cubic operator can be represented by a modal decomposition whose 
coefficients can easily be calculated numerically. Finally, let us note that, although we 
have illustrated our reasoning by considering the curvature perturbation ( only, it should 
be clear that similar techniques based on separability hold in theories with several degrees 
of freedom, like in multiple field inflationary models. 

2.3 Numerical considerations 

It turns out that most of the difficulties we encounter to numerically evaluate the Fourier 
space integrals needed to determine the modal decomposition coefficients can be under- 
stood by using the crudest approximation for the mode functions, namely 

Ofc(r) oc -L (1 + ikc.T)e-*°- r , C(r) oc -L^re"^ , (2.11) 
yk 6 yk 6 
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up to some slowly varying multiplicative factor function of both k and t, and where c s is 
the slowly- varying speed of sound of fluctuations. Required integrals of the type 

/*» ^ / 2fe-(U + U) \ k%(r)C k (r') dk (2.12) 

Jk m i n \ ^max #min / 

in (2.9) can be seem as sum of terms of the form 

^max 

/ k 2+n Ur)C k {r')dk (2.13) 

" ^min 

for n > 0. An obvious difficulty in numerically evaluating these integrals is the rapid 
oscillations coming from the phase factor e~ lkCsT in the mode functions. However, this 
can be circumvented by numerically solving the linear mode function equation not for 
Ck(f) itself, but rather for Cfc( r ) = Ck{T)e lkCsT , which displays no rapid oscillations. This 
results in the integrand in Eq. (2.13) being put in the form of the slowly varying function 
fc 2+n C/c( r )C( r/ ) multiplied by the explicitly known phase factor e L K c ^ T ') T> '- C 4 T ) T ) ? a type of 
integrand for which known efficient numerical routines can be used. In practice of course, 
the above slowly-varying function of k is only known as a result of evaluating the mode 
functions at sample /c-points and constructing an interpolating function out of them. This 
results in one last difficulty: for n — 0, k 2 + n (k(T)Q(r f ) contain terms behaving like 1/k 
which do not admit a rapidly convergent polynomial representation. We overcome this 
by simply considering the well behaved function /c 3 C/ c (r)C^(r / ) (instead of k 2 (k(T)Q(r f ) in 
Eq. (2.12)), and dividing by the appropriate factor of k at the end of the calculations. On 
the other hand, note that no such trick is required for evaluating the integrals of the form 
£r Pi fc2 Cfc( r )C'( r ') dk because of the additional k 2 factor in C k (r) in Eq. (2.11). 

Eventually, we observed that a slightly different expansion leads to improved accuracy, 
which consists in considering the shapes as function of the yi = log(^)(z = 1, 2, 3), rather 
than the k^s themselves, and therefore in decomposing them similarly as described above, 
but with the scalar product 

F log (S, S') = [ S( yi ,y 2 , ys)S'( yi , y 2 , y 3 ) d Vl dy 2 dy 3 . (2.14) 

J [log(/c min ),log(A; max I 

In the following, we denote results obtained by this method as resulting from the log- 
method. 

2.4 Projecting onto the interacting vacuum and the ie prescription 

When directly computing the value of a non-Gaussian shape in a specific momentum 
configuration, using e.g. Eq. (2.4), the ie prescription is needed to project on to the inter- 
acting vacuum. However, the quantities we are calculating in our approach are different: 
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these are the coefficients that enter into the shape's modal decomposition (2.5), like a n in 
Eq. (2.9). Here, the time-integrand comes as the result of integrating the mode functions in 
Fourier-space, which averages out the rapid oscillations of the mode functions deep inside 
the horizon, with one important consequence: evaluating Eq. (2.9) for e > and taking 
the limit e 0, following the ie prescription, is equivalent to just evaluating it at e = 
straight away. Or, to put it a different way: the process of first integrating over Fourier 
space before performing the required time-integral is equivalent to the ie prescription for 
us. 

Hence, while integrals like the one in Eq. (2.4) do not converge, as their lower bound 
7~ini goes to — oc, without an appropriate rotation of the time integration contour into the 
imaginary plane, this is not the case in our method. Actually, by using the asymptotic 
form of the mode functions (2.11) deep inside the horizon — which only results from our 
choice of the Bunch-Davies vacuum — and the behavior of the scale factor a{r) oc 1/r, 
one can easily determine the early-time behavior of the required time-integrands, like the 
one in Eq. (2.9). One thus finds that the coefficients of the modal decomposition converge, 
as Ti n i goes to — oc, like J°° e lx for the operators in J dta(t) ( 3 and J dta 3 (t) (( 2 (what- 
ever their explicit gradient structures), and like J°° ^e lx for the operators in J dta 3 (t) ( 3 
and J dta(t)(( 2 . In practice, these behaviors, which we have verified in our numerical 
calculations, imply that requisite time-integrations can be performed only starting from a 
few e-folds before the largest scale of interest crosses the horizon. 

3 Results 

In this section, we demonstrate that the method described above works as intended. We 
first describe the details of our set-up, in particular the inflationary backgrounds and the 
cubic operators we have considered, before showing explicit results in various situations. 

3.1 Set-up 

Test scenarios. — We consider two kinds of inflationary scenarios: the one of a canonical 
scalar field with an m 2 (j) 2 potential, which we refer to as the chaotic scenario, and the 
one corresponding to a Dirac-Born-Infeld (DBI) action in an AdS geometry [28] . In more 
detail, the two corresponding actions read (we omit the Einstein-Hilbert action which is 
common to the two set-up, and use units in which M p = 1) 

(3.1) 

l -m 2 <fA with m = A (3.2) 
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Our aim is not to consider inflationary scenarios that match precisely observational data, 
for example in terms of the amplitude and spectral index of the scalar fluctuations they 
generate or their level of non-Gaussianities, but rather to consider backgrounds in which 
all the slowly-varying approximations that are usually made in analytical treatments are 
satisfied, so that our numerical procedures can be tested against reliable theoretical pre- 
dictions. 

Unless otherwise stated, we take m = 10 -5 in both scenarios and A = 10 16 . We also 
use initial conditions such that background attractor solutions are reached rapidly, namely 

= — a/ §m and & = 10 3 in the chaotic case and <^ = — ( ^2 + ^ J and & = 2 in 
the DBI one. These parameters result in almost de-Sitter inflationary backgrounds in 
which the important background quantities e = —H/H 2 and c 2 = 1 — f((j))(f) 2 take al- 
most constant very small values during more than 60 efolds of expansion: e c h ao tic ~ 10 -6 , 
^dbi ~ 10 -3 and c s ~ 1CT 3 . 

As for the linear fluctuations about these inflationary backgrounds, the curvature 
perturbation obeys the simple second-order equation 

C + 2-C + c^ 2 a = 0, z=^, (3.3) 
z c s 

where it is understood that c 2 = 1 in the chaotic case, and we impose initial conditions 
such that the Bunch-Davies behavior holds inside the horizon: 



Ck = iH J ^-re- ikCsT for - kc s r > 1 . (3.4) 
V 4e/c 

Eventually, we consider a range of scales such that k ^.^ si = 10 3 — all studied scales are 
therefore well inside the horizon initially 1 — and with fc max //c min equals to 10 2 (similar 
results are obtained with fc max //c min = 10 3 ). 

Cubic operators. — Using the background and the linear fluctuations just described, we 
would like to test our method for numerically calculating the non-Gaussian shapes gen- 
erated by various cubic operators. As we have explained in section (2), we factor out 
any explicitly known fc-dependence coming from gradient terms in cubic operators. We 
therefore have to apply our method only to four kinds of cubic combinations of £, namely 
C 3 > C 2 C CC 2 an d C 3 > from which one can then reconstruct the shapes of any specific opera- 
tor. In the context of the most general second-order scalar-tensor theory, also known as 
Horndeski theory [29-31], it has been shown in Ref. [32] that only 5 independent cubic 



lr This is in agreement with our discussion about early-time convergence in subsection (2.4). 
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interactions can possibly be generated, i.e. that the third-order action can effectively be 
put in the form 



5 ( 3) = J dtd 3 xa 3 { 9l ( 3 + g 2 (( 2 + g 3 (^ + g^d* (d~ 2 c) +g 5 d\ (^" 2 c) 2 },(3.5) 

where the ^'s are time-dependent coupling constants, and with other interactions, like 
for example J dtd 3 xa((d() 2 , being redundant, in the sense that they can be expressed in 
terms of these 5 by using the linear equation of motion (3.3). This fact provides the op- 
portunity of interesting consistency checks, and we have therefore additionally considered 
this latter operator. 

Ordering of the basis functions. — We have described in subsection (2.2) our choice for the 
basis functions Q n we are using to expand any shape in a separable form. Of course, 
in any practical application, we should truncate the infinite expansion (2.5) and use a 
specific ordering for the Q n 's. For that purpose, we use the 'slicing' ordering described 
in Ref. [13] and decompose functions with polynomials (of ki or log(fcj) depending on the 
chosen method) up to a specific total order iV max . For symmetric functions for example, 
the index n thus relates to the completely symmetric triplets {ni,n 2 ,n 3 } (with a specific 
choice of sub-ordering): 

-> 000 4 -►111 8^ 022 12 — ^ 113 

1 -> 001 5^012 9^013 13^023 

2^011 6 -> 003 10 -> 004 14 ^014 (3.6) 

3 -> 002 7 ^ 112 11 ^ 122 15 ^ 005 ■ ■ ■ , 



where the transitions between total polynomial order is underlined. 
3.2 Numerical results 

Bare operators. — To show the effectiveness of our method, we applied it to the six op- 
erators listed above, with constant gr^'s in Eq. (3.5), for the two different test scenarios 
and using iV max = 6 in the modal decomposition 2 . We then evaluated the correlations — 
as measured by F(S, S')/ (^/F{S ) S)F(S I ', S'^J with the scalar product F (2.7) restricted 
to the tetrahedral domain — between the corresponding theoretical and numerically cal- 
culated shapes. Our results are listed in table (1). The various cosines are impressively 
close to 1, demonstrating the ability of our method to reconstruct non-Gaussian shapes. 

2 We stress that this is only methodological: for instance, not all 6 operators are present in the chaotic 
or DBI third-order action. 
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operator 


correlation (m 2 (p 2 ) 


correlation (DBI) 


c 3 


0.9992 


0.9994 


C(dC) 2 


0.99997 


0.99995 


cc 2 


0.999994 


0.999990 




0.999998 


0.999995 


^<r(0 -a 


0.99998 


0.99997 


dX&d-H? 


0.999990 


0.99998 



Table 1. Correlations between theoretical and numerically calculated shapes generated by var- 
ious (bare) operators on our two test background inflationary scenarios. 



Additionally, its rapid convergence is displayed in table (2), where we indicate the correla- 
tions between theoretical and numerically calculated shapes for increasing 7V max . We give 
the results for two representative operators of the equilateral and local type respectively, 
namely ( 3 and ((d() 2 , and for the test chaotic inflationary background (similar results 
are obtained in the DBI case). Eventually, let us note that it takes about 1.5 minutes 
on a laptop to calculate all possible bispectrum shapes (remember that they can be all 
deduced from the 4 cubic combinations of ( and () with 7V max = 2, and 4 minutes with 
^ m ax = 6. While this is sufficient for our purpose, we expect that a significant speed-up 
can be obtained after optimization of our code. 



operator 


iVmax = 


iVmax = 1 


N — 2 

1 "max " 


N max — 3 


N — 4 

1 "max ^ 


A^max — 5 


Aniax = 6 


c 3 


0.98 


0.97 


0.994 


0.998 


0.9990 


0.9993 


0.9994 




0.90 


0.997 


0.9997 


0.99994 


0.999990 


0.999996 


0.999998 



Table 2. Correlations between theoretical and numerically calculated shapes generated by the 
two bare operators £ 3 and C(d() 2 on the test chaotic background, for varying iV max - 

Two non-trivial checks. — In single-field inflation, there exists a well known consistency 
relation between the amplitude of the primordial bispectrum in the squeezed limit and the 
scalar spectral index [24, 25], namely that 

lim <C kl Ck 2 Ck 3 > = -M^H^Xn^) - l)P c {h)Pd^) (3-7) 

1 

where n 8 {k) is defined as 
with 

(C kl Ck 2 ) = (27r) 3 «5( 3 )(k 1 + k 2 )P c (*i) . (3.9) 
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To demonstrate the accuracy of our method, we calculated the bispectrum generated in 
various chaotic inflationary scenarios, which, at leading-order in a slow- varying approxi- 
mation, simply corresponds to setting g 2 (t) = e(e — rj) and gs(t) = e(e + rj) in Eq. (3.5), 
where rj = e/{He). We then evaluated the reconstructed bispectra on squeezed triangles 
such that ks = O.Olfci = 0.01/c 2 and we numerically calculated the corresponding n s {k\). 
We found that in all cases the single field consistency relation was verified to better than 
the % level. For instance, for (pi = 16, we find that n s (/c max ) — 1 = —0.0379, to be compared 
with —0.0377 deduced from the numerically calculated bispectrum, so a 0.5% discrepancy. 



As we mentioned above, another non-trivial check of our numerical method involves 
the redundancy of some cubic operators [32] . At leading order in the slow- varying approx- 
imation, the DBI third-order action reads 



S { 



(3),DBI 



dt d 6 x a e [ — — 1 



C(<9C) 5 



(3.10) 



Surprisingly enough at first sight, the two cubic operators present in this action generate 
non-Gaussian shapes that are strongly correlated with the local ansatz [34], whereas the 
DBI scenario is known to generate a bispectrum of the equilateral type [35]. A way of 
understanding this is to note that, upon using the linear mode function equation (3.3), an 
equivalent DBI third-order action reads, at leading-order in the slow- varying approxima- 
tion 3 [32]: 



(3),DBI 



~ / dt d 3 x a 3 



H 



1 



^c 3 + 



(3.11) 



This latter form of the third-order action makes it manifest that the total DBI bispectrum 
is of equilateral type, as this is the case for the individual bispectra generated by the 
operators in ( 3 and ((d() 2 . 

The fact that the local-type behaviors of the two operators (( 2 and ((d() 2 should 
conspire to cancel in (3.10), resulting in a total equilateral-type shape, provides a non- 
trivial test of our method. We have verified it, obtaining a correlation of 0.993 between 
the bispectra calculated using Eq. (3.10) and Eq. (3.11) in our test DBI scenario. This is 
summarized graphically in Fig. (1), where shapes are represented through several density 
contours. 



Other examples. — For illustrative purposes, we display below examples of theoretical bis- 
pectrum shapes (left) and their numerically calculated counterparts (right) in our test DBI 

3 Note that the first form (3.10) of naturally results from calculating it in the uniform inflaton 
gauge, while the form (3.11) naturally appears when working in the spatially flat gauge. 
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Figure 1. A graphical representation of a non-trivial check of our method. From left to right: 
the shapes generated in our test DBI scenario by, respectively, the (( 2 and C(^C) 2 operators in 
Eq. (3.10), the resulting total shape, and the total shape calculated from Eq. (3.11). Each shape 
are represented through several density contours. The third one shows that the cancellations of 
the local-type behaviors of the first two shapes is not perfect in the squeezed limits (in blue). 
The overlap with the fourth shape is nonetheless remarkable, with a correlation of 0.993 between 
the two. 

scenario. We used the log-method (c.f. subsection (2.3)) with iV max = 8, together with 
the coupling constants of Eqs. (3.10)-(3.11). The agreement is impressive. 

Theoretical bispectra Numerically calculated bispectra 
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Note that because the slow- varying approximation is very accurate on our test DBI 
background, the numerically calculated shapes displayed above are almost exactly scale 
invariant. For comparison, we show in Fig. (2) the shapes generated by the two operators 
C 3 and C(d() 2 in the third-order action Eq. (3.11) under the same conditions as previously, 
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but with A = 10 12 in the action (3.2). This modification results in a much more rapidly 
varying background, and hence to strongly scale-dependent equilateral non-Gaussianities, 
whose amplitude in the equilateral limit grows by a factor of about 30 between the largest 
and the smallest scales. This exemplifies the effectiveness of our approach for studying 
primordial non-Gaussianities generated from non-trivial inflationary dynamics. 




Figure 2. Shapes generated by the two operators £ 3 (left) and C(d£) 2 (right) in the third-order 
action Eq. (3.11), on the DBI background with A = 10 12 (see the main text). The rapidly varying 
background leads to strongly scale-dependent equilateral-type bispectra. 

4 Conclusions 

In this paper we have proposed a new method to numerically calculate higher-order cor- 
relation functions of primordial fluctuations generated from early-universe scenarios. It 
is based on the realization that the tree-level In-In formalism is intrinsically separable, 
which enables one to use modal techniques to efficiently represent any non-Gaussian pro- 
file. We have proved the feasibility and the accuracy of our method by applying it to 
simple single-field inflationary scenarios in which analytical results are available, finding 
impressive agreements and rapid convergence with respect to the number of modes used 
in the decomposition. We also performed non-trivial consistency checks, including the 
verification of the single field consistency relation, and the cancellations of two local-type 
behaviors to obtain an equilateral-type shape, as relevant in the DBI inflationary model. 
Additionally, we found our code to be quick, although a significant speed-up can be ex- 
pected after proper optimization. 
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To us, the main advantages of our method are that: (i) the ie prescription is auto- 
matically taken into account, preventing the need for ad-hoc tricks to implement it. (ii) 
We are able to calculate correlation functions at once over the whole domain of interest, 
and not only for specific fc-space configuration, (iii) Non-Gaussian shapes are directly 
expressed in a separable form, and hence are ready to be efficiently constrained by data. 

While we have mainly applied our method to simple early-universe models, thus 
providing a proof of its concept, its usefulness ultimately lies in the fact that it provides 
an efficient way of determining the non-Gaussian properties of cosmological fluctuations in 
non-trivial models in which analytical techniques are known to be insufficient. The result 
of such investigations will be presented elsewhere [36, 37]. 
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